Shibboleth model 1 analyses

#take averages

cost10.B0.avg <- trait.avg(s1.new.B0.cost10[[1]])
cost10.B.1.avg <- trait.avg(s1.new.B.1.cost10[[1]])
cost10.B.2.avg <- trait.avg(s1.new.B.2.cost10[[1]])
cost10.B.3.avg <- trait.avg(s1.new.B.3.cost10[[1]])
cost10.B.4.avg <- trait.avg(s1.new.B.4.cost10[[1]])
cost10.B.5.avg <- trait.avg(s1.new.B.5.cost10[[1]])
cost10.B.6.avg <- trait.avg(s1.new.B.6.cost10[[1]])
cost10.B.7.avg <- trait.avg(s1.new.B.7.cost10[[1]])
cost10.B.8.avg <- trait.avg(s1.new.B.8.cost10[[1]])
cost10.B.9.avg <- trait.avg(s1.new.B.9.cost10[[1]])
cost10.B1.avg <- trait.avg(s1.new.B1.cost10[[1]])

#group 

cost10.averages <- rbind(cost10.B0.avg, cost10.B.1.avg, cost10.B.2.avg, cost10.B.3.avg, cost10.B.4.avg, cost10.B.5.avg, cost10.B.6.avg, cost10.B.7.avg, cost10.B.8.avg, cost10.B.9.avg, cost10.B1.avg)

cost10.averages$B <- NA
cost10.averages$B <- rep(seq(0,1,0.1), each=2000)

library(ggplot2)
ggplot(data=subset(cost10.averages, generation==200), aes(x=B,y=Freq, color=Var2))+
geom_line()+
scale_color_manual(values=unique(cost10.averages$Var2))+labs(y="frequency")+theme_bw()+theme(legend.position="none")

cost10.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost10[[2]], s1.new.B.1.cost10[[2]], s1.new.B.2.cost10[[2]], s1.new.B.3.cost10[[2]], s1.new.B.4.cost10[[2]], s1.new.B.5.cost10[[2]], s1.new.B.6.cost10[[2]], s1.new.B.7.cost10[[2]], s1.new.B.8.cost10[[2]], s1.new.B.9.cost10[[2]], s1.new.B1.cost10[[2]])

cost10.B.drift.outcomes$B <- NA
cost10.B.drift.outcomes$B[1:5000] <- "drift"
cost10.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost10.B.drift.outcomes$B <- factor(cost10.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost10.B.drift.outcomes$green.outcome <- NA
cost10.B.drift.outcomes$green.outcome[cost10.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost10.B.drift.outcomes$green.outcome[cost10.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost10.model <- glm(green.outcome~B, data=cost10.B.drift.outcomes, family="binomial")
summary(cost10.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost10.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5764  -1.1603   0.8254   1.1946   1.5317  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.040397   0.028864  -1.400  0.16165    
B0           0.942299   0.115693   8.145 3.80e-16 ***
B0.1         0.759911   0.110382   6.884 5.80e-12 ***
B0.2         0.771284   0.110586   6.975 3.07e-12 ***
B0.3         0.488946   0.105539   4.633 3.61e-06 ***
B0.4         0.282783   0.105032   2.692  0.00710 ** 
B0.5         0.342678   0.105182   3.258  0.00112 ** 
B0.6        -0.004951   0.104470  -0.047  0.96220    
B0.7        -0.232382   0.106089  -2.190  0.02849 *  
B0.8        -0.308271   0.106288  -2.900  0.00373 ** 
B0.9        -0.762501   0.111433  -6.843 7.77e-12 ***
B1          -0.565129   0.105866  -5.338 9.39e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12759  on 9204  degrees of freedom
Residual deviance: 12437  on 9193  degrees of freedom
  (1295 observations deleted due to missingness)
AIC: 12461

Number of Fisher Scoring iterations: 4

plot_model(cost10.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

#determining likelihood of any trait going into fixation by generation 200

cost10.B.drift.outcomes$any.outcome <- NA
cost10.B.drift.outcomes$any.outcome[is.na(cost10.B.drift.outcomes$e)] <- FALSE
cost10.B.drift.outcomes$any.outcome[!is.na(cost10.B.drift.outcomes$e)] <- TRUE

#create model and plot

cost10.model.b <- glm(any.outcome~B, data=cost10.B.drift.outcomes, family="binomial")
summary(cost10.model.b)

Call:
glm(formula = any.outcome ~ B, family = "binomial", data = cost10.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5432   0.2835   0.2835   0.6681   0.7122  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.19379    0.07269   43.94   <2e-16 ***
B0          -1.95129    0.12958  -15.06   <2e-16 ***
B0.1        -1.80750    0.13336  -13.55   <2e-16 ***
B0.2        -1.80750    0.13336  -13.55   <2e-16 ***
B0.3        -1.70431    0.13640  -12.49   <2e-16 ***
B0.4        -1.83231    0.13267  -13.81   <2e-16 ***
B0.5        -1.80750    0.13336  -13.55   <2e-16 ***
B0.6        -1.84458    0.13233  -13.94   <2e-16 ***
B0.7        -1.91643    0.13045  -14.69   <2e-16 ***
B0.8        -1.88088    0.13136  -14.32   <2e-16 ***
B0.9        -1.75673    0.13482  -13.03   <2e-16 ***
B1          -1.50550    0.14309  -10.52   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7843.8  on 10499  degrees of freedom
Residual deviance: 7148.4  on 10488  degrees of freedom
AIC: 7172.4

Number of Fisher Scoring iterations: 6

plot_model(cost10.model.b, show.values = TRUE, value.offset=.3, title="Odds ratios for any outcome")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

#B.25

cost0.B.25.avg <- trait.avg(s1.new.B.25.cost0[[1]])
cost10.B.25.avg <- trait.avg(s1.new.B.25.cost10[[1]])
cost20.B.25.avg <- trait.avg(s1.new.B.25.cost20[[1]])
cost30.B.25.avg <- trait.avg(s1.new.B.25.cost30[[1]])
cost40.B.25.avg <- trait.avg(s1.new.B.25.cost40[[1]])
cost50.B.25.avg <- trait.avg(s1.new.B.25.cost50[[1]])
cost60.B.25.avg <- trait.avg(s1.new.B.25.cost60[[1]])
cost70.B.25.avg <- trait.avg(s1.new.B.25.cost70[[1]])
cost80.B.25.avg <- trait.avg(s1.new.B.25.cost80[[1]])
cost90.B.25.avg <- trait.avg(s1.new.B.25.cost90[[1]])
cost100.B.25.avg <- trait.avg(s1.new.B.25.cost100[[1]])

#group 

B.25.averages <- rbind(cost0.B.25.avg, cost10.B.25.avg, cost20.B.25.avg, cost30.B.25.avg, cost40.B.25.avg, cost50.B.25.avg, cost60.B.25.avg, cost70.B.25.avg, cost80.B.25.avg, cost90.B.25.avg, cost100.B.25.avg)

B.25.averages$cost <- NA
B.25.averages$cost <- rep(seq(0,100,10), each=2000)

library(ggplot2)
ggplot(data=subset(B.25.averages, generation==200), aes(x=cost,y=Freq, color=Var2))+
geom_line()+
scale_color_manual(values=unique(B.25.averages$Var2))+labs(y="frequency")+theme_bw()+theme(legend.position="none")

B.25.cost.drift.outcomes <- rbind(drift5000[[2]], s1.new.B.25.cost0[[2]], s1.new.B.25.cost10[[2]], s1.new.B.25.cost20[[2]], s1.new.B.25.cost30[[2]], s1.new.B.25.cost40[[2]], s1.new.B.25.cost50[[2]], s1.new.B.25.cost60[[2]], s1.new.B.25.cost70[[2]], s1.new.B.25.cost80[[2]], s1.new.B.25.cost90[[2]], s1.new.B.25.cost100[[2]])

B.25.cost.drift.outcomes$cost <- NA
B.25.cost.drift.outcomes$cost[1:5000] <- "drift"
B.25.cost.drift.outcomes$cost[5001:10500] <- rep(seq(0,100,10), each=500)
B.25.cost.drift.outcomes$cost <- factor(B.25.cost.drift.outcomes$cost, levels=c("drift","0","10","20","30","40","50","60","70","80","90","100"))

library(colorRamps)
B.25.cost.drift.outcomes$green.outcome <- NA
B.25.cost.drift.outcomes$green.outcome[B.25.cost.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
B.25.cost.drift.outcomes$green.outcome[B.25.cost.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

summary(B.25.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B.25.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5898  -1.1603   0.8383   1.1946   1.1946  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
cost0        0.33532    0.10477   3.201  0.00137 ** 
cost10       0.75274    0.11130   6.763 1.35e-11 ***
cost20       0.87331    0.11309   7.722 1.14e-14 ***
cost30       0.60587    0.10758   5.632 1.78e-08 ***
cost40       0.97196    0.11290   8.609  < 2e-16 ***
cost50       0.68635    0.10790   6.361 2.01e-10 ***
cost60       0.90539    0.11258   8.042 8.84e-16 ***
cost70       0.66807    0.10786   6.194 5.88e-10 ***
cost80       0.93899    0.11479   8.180 2.84e-16 ***
cost90       0.79641    0.10978   7.255 4.02e-13 ***
cost100      0.81483    0.11232   7.255 4.03e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12578  on 9230  degrees of freedom
Residual deviance: 12235  on 9219  degrees of freedom
  (1269 observations deleted due to missingness)
AIC: 12259

Number of Fisher Scoring iterations: 4

plot_model(B.25.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

B.25.model.b <- glm(any.outcome~cost, data=B.25.cost.drift.outcomes, family="binomial")
summary(B.25.model.b)

Call:
glm(formula = any.outcome ~ cost, family = "binomial", data = B.25.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5432   0.2835   0.2835   0.6492   0.6976  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.19379    0.07269   43.94   <2e-16 ***
cost0       -1.76957    0.13444  -13.16   <2e-16 ***
cost10      -1.90466    0.13075  -14.57   <2e-16 ***
cost20      -1.85677    0.13200  -14.07   <2e-16 ***
cost30      -1.76957    0.13444  -13.16   <2e-16 ***
cost40      -1.62227    0.13902  -11.67   <2e-16 ***
cost50      -1.67744    0.13724  -12.22   <2e-16 ***
cost60      -1.74378    0.13520  -12.90   <2e-16 ***
cost70      -1.70431    0.13640  -12.49   <2e-16 ***
cost80      -1.88088    0.13136  -14.32   <2e-16 ***
cost90      -1.67744    0.13724  -12.22   <2e-16 ***
cost100     -1.89281    0.13105  -14.44   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7741.2  on 10499  degrees of freedom
Residual deviance: 7078.6  on 10488  degrees of freedom
AIC: 7102.6

Number of Fisher Scoring iterations: 6

plot_model(B.25.model.b, show.values = TRUE, value.offset=.3, title="Odds ratios for any outcome")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

#B.6

cost0.B.6.avg <- trait.avg(s1.new.B.6.cost0[[1]])
cost10.B.6.avg <- trait.avg(s1.new.B.6.cost10[[1]])

#group 

B.6.averages <- rbind(cost0.B.6.avg, cost10.B.6.avg, cost20.B.6.avg, cost30.B.6.avg, cost40.B.6.avg, cost50.B.6.avg, cost60.B.6.avg, cost70.B.6.avg, cost80.B.6.avg, cost90.B.6.avg, cost100.B.6.avg)

B.6.averages$cost <- NA
B.6.averages$cost <- rep(seq(0,100,10), each=2000)

library(ggplot2)
ggplot(data=subset(B.6.averages, generation==200), aes(x=cost,y=Freq, color=Var2))+
geom_line()+
scale_color_manual(values=unique(B.6.averages$Var2))+labs(y="frequency")+theme_bw()+theme(legend.position="none")

B.6.cost.drift.outcomes <- rbind(drift5000[[2]], s1.new.B.6.cost0[[2]], s1.new.B.6.cost10[[2]], s1.new.B.6.cost20[[2]], s1.new.B.6.cost30[[2]], s1.new.B.6.cost40[[2]], s1.new.B.6.cost50[[2]], s1.new.B.6.cost60[[2]], s1.new.B.6.cost70[[2]], s1.new.B.6.cost80[[2]], s1.new.B.6.cost90[[2]], s1.new.B.6.cost100[[2]])

B.6.cost.drift.outcomes$cost <- NA
B.6.cost.drift.outcomes$cost[1:5000] <- "drift"
B.6.cost.drift.outcomes$cost[5001:10500] <- rep(seq(0,100,10), each=500)
B.6.cost.drift.outcomes$cost <- factor(B.6.cost.drift.outcomes$cost, levels=c("drift","0","10","20","30","40","50","60","70","80","90","100"))

library(colorRamps)
B.6.cost.drift.outcomes$green.outcome <- NA
B.6.cost.drift.outcomes$green.outcome[B.6.cost.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
B.6.cost.drift.outcomes$green.outcome[B.6.cost.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
B.6.model <- glm(green.outcome~cost, data=B.6.cost.drift.outcomes, family="binomial")
summary(B.6.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B.6.cost.drift.outcomes)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.244  -1.160  -1.069   1.195   1.289  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.040397   0.028864  -1.400   0.1617  
cost0       -0.219280   0.105504  -2.078   0.0377 *
cost10      -0.048771   0.103701  -0.470   0.6381  
cost20       0.196099   0.104494   1.877   0.0606 .
cost30      -0.008874   0.103399  -0.086   0.9316  
cost40      -0.061213   0.104936  -0.583   0.5597  
cost50      -0.044020   0.103809  -0.424   0.6715  
cost60       0.135231   0.104070   1.299   0.1938  
cost70      -0.057242   0.103017  -0.556   0.5784  
cost80       0.050153   0.102905   0.487   0.6260  
cost90      -0.035976   0.105005  -0.343   0.7319  
cost100      0.116000   0.104514   1.110   0.2670  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12771  on 9214  degrees of freedom
Residual deviance: 12758  on 9203  degrees of freedom
  (1285 observations deleted due to missingness)
AIC: 12782

Number of Fisher Scoring iterations: 3

plot_model(B.6.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

#determining likelihood of any trait going into fixation by generation 200

B.6.cost.drift.outcomes$any.outcome <- NA
B.6.cost.drift.outcomes$any.outcome[is.na(B.6.cost.drift.outcomes$e)] <- FALSE
B.6.cost.drift.outcomes$any.outcome[!is.na(B.6.cost.drift.outcomes$e)] <- TRUE

#create model and plot

B.6.model.b <- glm(any.outcome~cost, data=B.6.cost.drift.outcomes, family="binomial")
summary(B.6.model.b)

Call:
glm(formula = any.outcome ~ cost, family = "binomial", data = B.6.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5432   0.2835   0.2835   0.6568   0.6940  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.19379    0.07269   43.94   <2e-16 ***
cost0       -1.86887    0.13168  -14.19   <2e-16 ***
cost10      -1.75673    0.13482  -13.03   <2e-16 ***
cost20      -1.81995    0.13301  -13.68   <2e-16 ***
cost30      -1.73073    0.13559  -12.76   <2e-16 ***
cost40      -1.88088    0.13136  -14.32   <2e-16 ***
cost50      -1.76957    0.13444  -13.16   <2e-16 ***
cost60      -1.79495    0.13371  -13.42   <2e-16 ***
cost70      -1.67744    0.13724  -12.22   <2e-16 ***
cost80      -1.67744    0.13724  -12.22   <2e-16 ***
cost90      -1.89281    0.13105  -14.44   <2e-16 ***
cost100     -1.84458    0.13233  -13.94   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7804.5  on 10499  degrees of freedom
Residual deviance: 7126.7  on 10488  degrees of freedom
AIC: 7150.7

Number of Fisher Scoring iterations: 6

plot_model(B.6.model.b, show.values = TRUE, value.offset=.3, title="Odds ratios for any outcome")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

#extract odds ratios and confidence intervals from 5 models — varying B and holding cost at 10

B.model1.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(cost10.model$coefficients)[-1], low = exp(confint(cost10.model))[-1,1], high = exp(confint(cost10.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model2.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(s2.cost10.model$coefficients)[-1], low = exp(confint(s2.cost10.model))[-1,1], high = exp(confint(s2.cost10.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model3.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(s3.cost.model$coefficients)[-1], low = exp(confint(s3.cost.model))[-1,1], high = exp(confint(s3.cost.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model4.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(s4.cost.model$coefficients)[-1], low = exp(confint(s4.cost.model))[-1,1], high = exp(confint(s4.cost.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model5.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(s5.cost.model$coefficients)[-1], low = exp(confint(s5.cost.model))[-1,1], high = exp(confint(s5.cost.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model.outputs <- rbind(B.model1.outputs, B.model2.outputs, B.model3.outputs, B.model4.outputs, B.model5.outputs)

ggplot(data = B.model.outputs, aes(x = as.numeric(as.character(B)), y = OR, color = model))+
geom_line()+
geom_point(size = 1)+
geom_ribbon(aes(y = OR, ymin = low, ymax = high), alpha = 0.2)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
ylim(0, 111)+
ylab("Odds ratios for a green outcome")+
xlab("Value of B")+
scale_x_continuous(breaks = seq(0, 1, 0.1)) 

ggplot(data = B.model.outputs, aes(x = as.numeric(as.character(B)), y = OR, color = model))+
geom_line()+
geom_point(size = 1)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
ylim(0, 111)+
ylab("Odds ratios for a green outcome")+
xlab("Value of B")+
scale_x_continuous(breaks = seq(0, 1, 0.1)) 

ggplot(data = B.model.outputs[-c(which(B.model.outputs$B == 0), which(B.model.outputs$B == 0.1), which(B.model.outputs$B == 0.2), which(B.model.outputs$B == 0.3)),], aes(x = as.numeric(as.character(B)), y = OR, color = model, fill = model))+
geom_line()+
geom_point(size = 1.2)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
geom_ribbon(aes(y = OR, ymin = low, ymax = high), alpha = 0.2)+
ylab("Odds ratios for a green outcome")+
xlab("Value of B")+
scale_x_continuous(breaks = seq(0.4, 1, 0.1)) 

ggplot(data = B.model.outputs[c(which(B.model.outputs$B == 0), which(B.model.outputs$B == 0.1), which(B.model.outputs$B == 0.2), which(B.model.outputs$B == 0.3)),], aes(x = as.numeric(as.character(B)), y = OR, color = model, fill = model))+
geom_line()+
geom_point(size = 1.2)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
geom_ribbon(aes(y = OR, ymin = low, ymax = high), alpha = 0.2)+
ylab("Odds ratios for a green outcome")+
xlab("Value of B")+
ylim(0,10)+
scale_x_continuous(breaks = seq(0, 0.4, 0.1)) 

#extract odds ratios and confidence intervals from 5 models — varying cost and holding B at 0.6 (0.7 and 0.8 in models 2 and 5, respectively)

cost.model1.outputs <- data.frame(model = rep(NA, 11), B = rep(0.6,11), cost = seq(0,100,10), OR = exp(B.6.model$coefficients)[-1], low = exp(confint(B.6.model))[-1,1], high = exp(confint(B.6.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model2.outputs <- data.frame(model = rep(NA, 11), B = rep(0.7,11), cost = seq(0,100,10), OR = exp(s2.B.7.model$coefficients)[-1], low = exp(confint(s2.B.7.model))[-1,1], high = exp(confint(s2.B.7.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model3.outputs <- data.frame(model = rep(NA, 11), B = rep(0.6,11), cost = seq(0,100,10), OR = exp(s3.B.6.model$coefficients)[-1], low = exp(confint(s3.B.6.model))[-1,1], high = exp(confint(s3.B.6.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model4.outputs <- data.frame(model = rep(NA, 11), B = rep(0.6,11), cost = seq(0,100,10), OR = exp(s4.B.6.model$coefficients)[-1], low = exp(confint(s4.B.6.model))[-1,1], high = exp(confint(s4.B.6.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model5.outputs <- data.frame(model = rep(NA, 11), B = rep(0.8,11), cost = seq(0,100,10), OR = exp(s5.B.8.model$coefficients)[-1], low = exp(confint(s5.B.8.model))[-1,1], high = exp(confint(s5.B.8.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model.outputs <- rbind(cost.model1.outputs, cost.model2.outputs, cost.model3.outputs, cost.model4.outputs, cost.model5.outputs)

cost.model.outputs$model[1:11] <- "Model 1"
cost.model.outputs$model[12:22] <- "Model 2"
cost.model.outputs$model[23:33] <- "Model 3"
cost.model.outputs$model[34:44] <- "Model 4"
cost.model.outputs$model[45:55] <- "Model 5"
cost.model.outputs$model <- factor(B.model.outputs$model)

ggplot(data = cost.model.outputs, aes(x = as.numeric(as.character(cost)), y = OR, color = model))+
geom_line()+
geom_point(size = 1)+
geom_ribbon(aes(y = OR, ymin = low, ymax = high, fill=model), alpha = 0.2)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
ylim(0, 2.5)+
ylab("Odds ratios for a green outcome")+
xlab("Value of cost")+
scale_x_continuous(breaks = seq(0, 100, 10)) 

B0.cost.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost0[[2]], s1.new.B0.cost10[[2]], s1.new.B0.cost20[[2]], s1.new.B0.cost30[[2]], s1.new.B0.cost40[[2]], s1.new.B0.cost50[[2]], s1.new.B0.cost60[[2]], s1.new.B0.cost70[[2]], s1.new.B0.cost80[[2]], s1.new.B0.cost90[[2]], s1.new.B0.cost100[[2]])

B0.cost.drift.outcomes$cost <- NA
B0.cost.drift.outcomes$cost[1:5000] <- "drift"
B0.cost.drift.outcomes$cost[5001:10500] <- rep(seq(0,100,10), each=500)
B0.cost.drift.outcomes$cost <- factor(B0.cost.drift.outcomes$cost, levels=c("drift","0","10","20","30","40","50","60","70","80","90","100"))

library(colorRamps)
B0.cost.drift.outcomes$green.outcome <- NA
B0.cost.drift.outcomes$green.outcome[B0.cost.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
B0.cost.drift.outcomes$green.outcome[B0.cost.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
B0.model <- glm(green.outcome~cost, data=B0.cost.drift.outcomes, family="binomial")
summary(B0.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B0.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6515  -1.1603   0.8024   1.1946   1.1946  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400    0.162    
cost0        0.13055    0.09930   1.315    0.189    
cost10       0.96368    0.11441   8.423  < 2e-16 ***
cost20       0.97704    0.11765   8.305  < 2e-16 ***
cost30       0.89809    0.11141   8.061 7.54e-16 ***
cost40       1.00853    0.11615   8.683  < 2e-16 ***
cost50       1.10886    0.11889   9.327  < 2e-16 ***
cost60       0.87918    0.11389   7.720 1.17e-14 ***
cost70       0.92423    0.11273   8.198 2.44e-16 ***
cost80       0.93737    0.11435   8.198 2.45e-16 ***
cost90       1.02796    0.11488   8.948  < 2e-16 ***
cost100      0.86202    0.11155   7.728 1.10e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12521  on 9241  degrees of freedom
Residual deviance: 12046  on 9230  degrees of freedom
  (1258 observations deleted due to missingness)
AIC: 12070

Number of Fisher Scoring iterations: 4

plot_model(B0.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

B.75.cost.drift.outcomes <- rbind(drift5000[[2]], s1.new.B.75.cost0[[2]], s1.new.B.75.cost10[[2]], s1.new.B.75.cost20[[2]], s1.new.B.75.cost30[[2]], s1.new.B.75.cost40[[2]], s1.new.B.75.cost50[[2]], s1.new.B.75.cost60[[2]], s1.new.B.75.cost70[[2]], s1.new.B.75.cost80[[2]], s1.new.B.75.cost90[[2]], s1.new.B.75.cost100[[2]])

B.75.cost.drift.outcomes$cost <- NA
B.75.cost.drift.outcomes$cost[1:5000] <- "drift"
B.75.cost.drift.outcomes$cost[5001:10500] <- rep(seq(0,100,10), each=500)
B.75.cost.drift.outcomes$cost <- factor(B.75.cost.drift.outcomes$cost, levels=c("drift","0","10","20","30","40","50","60","70","80","90","100"))

library(colorRamps)
B.75.cost.drift.outcomes$green.outcome <- NA
B.75.cost.drift.outcomes$green.outcome[B.75.cost.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
B.75.cost.drift.outcomes$green.outcome[B.75.cost.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
B.75.model <- glm(green.outcome~cost, data=B.75.cost.drift.outcomes, family="binomial")
summary(B.75.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B.75.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1795  -1.1603  -0.9685   1.1946   1.4357  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
cost0       -0.54890    0.10682  -5.139 2.77e-07 ***
cost10      -0.19994    0.10564  -1.893  0.05840 .  
cost20      -0.30294    0.10550  -2.871  0.00409 ** 
cost30      -0.32780    0.10559  -3.104  0.00191 ** 
cost40      -0.32305    0.10629  -3.039  0.00237 ** 
cost50      -0.43204    0.10665  -4.051 5.10e-05 ***
cost60      -0.32733    0.10620  -3.082  0.00206 ** 
cost70      -0.32968    0.10586  -3.114  0.00184 ** 
cost80      -0.45318    0.10690  -4.239 2.24e-05 ***
cost90      -0.47311    0.10752  -4.400 1.08e-05 ***
cost100      0.04538    0.10396   0.437  0.66244    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12664  on 9199  degrees of freedom
Residual deviance: 12577  on 9188  degrees of freedom
  (1300 observations deleted due to missingness)
AIC: 12601

Number of Fisher Scoring iterations: 4

plot_model(B.75.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

B.75.model <- glm(green.outcome~cost, data=B.75.cost.drift.outcomes, family="binomial")
summary(B.75.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B.75.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1795  -1.1603  -0.9685   1.1946   1.4357  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
cost0       -0.54890    0.10682  -5.139 2.77e-07 ***
cost10      -0.19994    0.10564  -1.893  0.05840 .  
cost20      -0.30294    0.10550  -2.871  0.00409 ** 
cost30      -0.32780    0.10559  -3.104  0.00191 ** 
cost40      -0.32305    0.10629  -3.039  0.00237 ** 
cost50      -0.43204    0.10665  -4.051 5.10e-05 ***
cost60      -0.32733    0.10620  -3.082  0.00206 ** 
cost70      -0.32968    0.10586  -3.114  0.00184 ** 
cost80      -0.45318    0.10690  -4.239 2.24e-05 ***
cost90      -0.47311    0.10752  -4.400 1.08e-05 ***
cost100      0.04538    0.10396   0.437  0.66244    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12664  on 9199  degrees of freedom
Residual deviance: 12577  on 9188  degrees of freedom
  (1300 observations deleted due to missingness)
AIC: 12601

Number of Fisher Scoring iterations: 4

plot_model(B.75.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

B1.cost.drift.outcomes <- rbind(drift5000[[2]], s1.new.B1.cost0[[2]], s1.new.B1.cost10[[2]], s1.new.B1.cost20[[2]], s1.new.B1.cost30[[2]], s1.new.B1.cost40[[2]], s1.new.B1.cost50[[2]], s1.new.B1.cost60[[2]], s1.new.B1.cost70[[2]], s1.new.B1.cost80[[2]], s1.new.B1.cost90[[2]], s1.new.B1.cost100[[2]])

B1.cost.drift.outcomes$cost <- NA
B1.cost.drift.outcomes$cost[1:5000] <- "drift"
B1.cost.drift.outcomes$cost[5001:10500] <- rep(seq(0,100,10), each=500)
B1.cost.drift.outcomes$cost <- factor(B1.cost.drift.outcomes$cost, levels=c("drift","0","10","20","30","40","50","60","70","80","90","100"))

library(colorRamps)
B1.cost.drift.outcomes$green.outcome <- NA
B1.cost.drift.outcomes$green.outcome[B1.cost.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
B1.cost.drift.outcomes$green.outcome[B1.cost.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
B1.model <- glm(green.outcome~cost, data=B1.cost.drift.outcomes, family="binomial")
summary(B1.model)

Call:
glm(formula = green.outcome ~ cost, family = "binomial", data = B1.cost.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1603  -1.1603  -0.9225   1.1946   1.4791  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
cost0       -0.48750    0.10659  -4.574 4.80e-06 ***
cost10      -0.64553    0.10794  -5.981 2.22e-09 ***
cost20      -0.46136    0.10579  -4.361 1.29e-05 ***
cost30      -0.30791    0.10482  -2.938  0.00331 ** 
cost40      -0.53970    0.10766  -5.013 5.36e-07 ***
cost50      -0.51641    0.10542  -4.898 9.66e-07 ***
cost60      -0.59379    0.11187  -5.308 1.11e-07 ***
cost70      -0.48570    0.10503  -4.624 3.76e-06 ***
cost80      -0.32593    0.10533  -3.094  0.00197 ** 
cost90      -0.63061    0.10892  -5.790 7.06e-09 ***
cost100     -0.58401    0.10757  -5.429 5.67e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12688  on 9281  degrees of freedom
Residual deviance: 12532  on 9270  degrees of freedom
  (1218 observations deleted due to missingness)
AIC: 12556

Number of Fisher Scoring iterations: 4

plot_model(B1.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

#model 1 B model outputs

cost.model1.B0.outputs <- data.frame(model = rep(NA, 11), B = rep(0,11), cost = seq(0,100,10), OR = exp(B0.model$coefficients)[-1], low = exp(confint(B0.model))[-1,1], high = exp(confint(B0.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model1.B.25.outputs <- data.frame(model = rep(NA, 11), B = rep(0.25,11), cost = seq(0,100,10), OR = exp(B.25.model$coefficients)[-1], low = exp(confint(B.25.model))[-1,1], high = exp(confint(B.25.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model1.B.6.outputs <- data.frame(model = rep(NA, 11), B = rep(0.6,11), cost = seq(0,100,10), OR = exp(B.6.model$coefficients)[-1], low = exp(confint(B.6.model))[-1,1], high = exp(confint(B.6.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model1.B.75.outputs <- data.frame(model = rep(NA, 11), B = rep(0.75,11), cost = seq(0,100,10), OR = exp(B.75.model$coefficients)[-1], low = exp(confint(B.75.model))[-1,1], high = exp(confint(B.75.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model1.B1.outputs <- data.frame(model = rep(NA, 11), B = rep(1,11), cost = seq(0,100,10), OR = exp(B1.model$coefficients)[-1], low = exp(confint(B1.model))[-1,1], high = exp(confint(B1.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.outputs <- rbind(cost.model1.B0.outputs, cost.model1.B.25.outputs, cost.model1.B.6.outputs, cost.model1.B.75.outputs, cost.model1.B1.outputs)

B.model1.outputs$model[1:11] <- "B = 0"
B.model1.outputs$model[12:22] <- "B = .25"
B.model1.outputs$model[23:33] <- "B = .6"
B.model1.outputs$model[34:44] <- "B = .75"
B.model1.outputs$model[45:55] <- "B = 1"
B.model1.outputs$model <- factor(B.model1.outputs$model)

ggplot(data = B.model1.outputs, aes(x = as.numeric(as.character(cost)), y = OR, color = model))+
geom_line()+
geom_point(size = 1)+
geom_ribbon(aes(y = OR, ymin = low, ymax = high, fill=model), alpha = 0.2)+
scale_color_manual(values=c(blue2green(5)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
ylim(0, 4)+
ylab("Odds ratios for a green outcome")+
xlab("Value of cost")+
scale_x_continuous(breaks = seq(0, 100, 10))+
theme(legend.title=element_blank())

cost0.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost0[[2]], s1.new.B.1.cost0[[2]], s1.new.B.2.cost0[[2]], s1.new.B.3.cost0[[2]], s1.new.B.4.cost0[[2]], s1.new.B.5.cost0[[2]], s1.new.B.6.cost0[[2]], s1.new.B.7.cost0[[2]], s1.new.B.8.cost0[[2]], s1.new.B.9.cost0[[2]], s1.new.B1.cost0[[2]])

cost0.B.drift.outcomes$B <- NA
cost0.B.drift.outcomes$B[1:5000] <- "drift"
cost0.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost0.B.drift.outcomes$B <- factor(cost0.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost0.B.drift.outcomes$green.outcome <- NA
cost0.B.drift.outcomes$green.outcome[cost0.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost0.B.drift.outcomes$green.outcome[cost0.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost0.model <- glm(green.outcome~B, data=cost0.B.drift.outcomes, family="binomial")
summary(cost0.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost0.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3511  -1.1603  -0.9141   1.1946   1.4658  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
B0           0.13055    0.09930   1.315  0.18862    
B0.1         0.33374    0.10369   3.219  0.00129 ** 
B0.2         0.43993    0.10353   4.249 2.14e-05 ***
B0.3         0.12166    0.10489   1.160  0.24610    
B0.4        -0.20399    0.10444  -1.953  0.05080 .  
B0.5         0.14368    0.10338   1.390  0.16457    
B0.6        -0.21928    0.10551  -2.078  0.03767 *  
B0.7        -0.47170    0.10497  -4.494 7.00e-06 ***
B0.8        -0.61612    0.10844  -5.682 1.33e-08 ***
B0.9        -0.56194    0.10761  -5.222 1.77e-07 ***
B1          -0.48750    0.10659  -4.574 4.80e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12897  on 9318  degrees of freedom
Residual deviance: 12746  on 9307  degrees of freedom
  (1181 observations deleted due to missingness)
AIC: 12770

Number of Fisher Scoring iterations: 4

plot_model(cost0.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

cost25.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost25[[2]], s1.new.B.1.cost25[[2]], s1.new.B.2.cost25[[2]], s1.new.B.3.cost25[[2]], s1.new.B.4.cost25[[2]], s1.new.B.5.cost25[[2]], s1.new.B.6.cost25[[2]], s1.new.B.7.cost25[[2]], s1.new.B.8.cost25[[2]], s1.new.B.9.cost25[[2]], s1.new.B1.cost25[[2]])

cost25.B.drift.outcomes$B <- NA
cost25.B.drift.outcomes$B[1:5000] <- "drift"
cost25.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost25.B.drift.outcomes$B <- factor(cost25.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost25.B.drift.outcomes$green.outcome <- NA
cost25.B.drift.outcomes$green.outcome[cost25.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost25.B.drift.outcomes$green.outcome[cost25.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost25.model <- glm(green.outcome~B, data=cost25.B.drift.outcomes, family="binomial")
summary(cost25.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost25.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5981  -1.1603   0.8087   1.1946   1.3994  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400 0.161651    
B0           0.99035    0.11544   8.579  < 2e-16 ***
B0.1         0.90721    0.11435   7.933 2.13e-15 ***
B0.2         0.94091    0.11429   8.233  < 2e-16 ***
B0.3         0.74121    0.11110   6.672 2.53e-11 ***
B0.4         0.39114    0.10448   3.743 0.000181 ***
B0.5         0.23097    0.10452   2.210 0.027115 *  
B0.6         0.18457    0.10397   1.775 0.075875 .  
B0.7        -0.22959    0.10557  -2.175 0.029654 *  
B0.8        -0.24444    0.10387  -2.353 0.018603 *  
B0.9        -0.35891    0.10528  -3.409 0.000651 ***
B1          -0.46783    0.10599  -4.414 1.01e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12760  on 9215  degrees of freedom
Residual deviance: 12433  on 9204  degrees of freedom
  (1284 observations deleted due to missingness)
AIC: 12457

Number of Fisher Scoring iterations: 4

plot_model(cost25.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

cost50.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost50[[2]], s1.new.B.1.cost50[[2]], s1.new.B.2.cost50[[2]], s1.new.B.3.cost50[[2]], s1.new.B.4.cost50[[2]], s1.new.B.5.cost50[[2]], s1.new.B.6.cost50[[2]], s1.new.B.7.cost50[[2]], s1.new.B.8.cost50[[2]], s1.new.B.9.cost50[[2]], s1.new.B1.cost50[[2]])

cost50.B.drift.outcomes$B <- NA
cost50.B.drift.outcomes$B[1:5000] <- "drift"
cost50.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost50.B.drift.outcomes$B <- factor(cost50.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost50.B.drift.outcomes$green.outcome <- NA
cost50.B.drift.outcomes$green.outcome[cost50.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost50.B.drift.outcomes$green.outcome[cost50.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost50.model <- glm(green.outcome~B, data=cost50.B.drift.outcomes, family="binomial")
summary(cost50.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost50.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6515  -1.1603   0.7685   1.1946   1.4211  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400 0.161651    
B0           1.10886    0.11889   9.327  < 2e-16 ***
B0.1         0.88648    0.11377   7.792 6.61e-15 ***
B0.2         0.66113    0.11018   6.001 1.96e-09 ***
B0.3         0.49359    0.10834   4.556 5.21e-06 ***
B0.4         0.25569    0.10452   2.446 0.014432 *  
B0.5         0.36937    0.10407   3.549 0.000387 ***
B0.6        -0.04402    0.10381  -0.424 0.671529    
B0.7        -0.33909    0.10784  -3.144 0.001664 ** 
B0.8        -0.09447    0.10418  -0.907 0.364516    
B0.9        -0.34649    0.10552  -3.284 0.001025 ** 
B1          -0.51641    0.10542  -4.898 9.66e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12731  on 9188  degrees of freedom
Residual deviance: 12441  on 9177  degrees of freedom
  (1311 observations deleted due to missingness)
AIC: 12465

Number of Fisher Scoring iterations: 4

plot_model(cost50.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

cost75.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost75[[2]], s1.new.B.1.cost75[[2]], s1.new.B.2.cost75[[2]], s1.new.B.3.cost75[[2]], s1.new.B.4.cost75[[2]], s1.new.B.5.cost75[[2]], s1.new.B.6.cost75[[2]], s1.new.B.7.cost75[[2]], s1.new.B.8.cost75[[2]], s1.new.B.9.cost75[[2]], s1.new.B1.cost75[[2]])

cost75.B.drift.outcomes$B <- NA
cost75.B.drift.outcomes$B[1:5000] <- "drift"
cost75.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost75.B.drift.outcomes$B <- factor(cost75.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost75.B.drift.outcomes$green.outcome <- NA
cost75.B.drift.outcomes$green.outcome[cost75.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost75.B.drift.outcomes$green.outcome[cost75.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost75.model <- glm(green.outcome~B, data=cost75.B.drift.outcomes, family="binomial")
summary(cost75.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost75.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6636  -1.1603   0.7597   1.1946   1.4404  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400 0.161651    
B0           1.13560    0.12021   9.447  < 2e-16 ***
B0.1         0.99522    0.11488   8.663  < 2e-16 ***
B0.2         0.91727    0.11284   8.129 4.33e-16 ***
B0.3         0.69744    0.10771   6.475 9.48e-11 ***
B0.4         0.27796    0.10404   2.672 0.007549 ** 
B0.5         0.20960    0.10553   1.986 0.047019 *  
B0.6         0.04040    0.10384   0.389 0.697262    
B0.7        -0.02180    0.10586  -0.206 0.836846    
B0.8        -0.30617    0.10457  -2.928 0.003412 ** 
B0.9        -0.55949    0.10697  -5.231 1.69e-07 ***
B1          -0.41280    0.10834  -3.810 0.000139 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12748  on 9205  degrees of freedom
Residual deviance: 12387  on 9194  degrees of freedom
  (1294 observations deleted due to missingness)
AIC: 12411

Number of Fisher Scoring iterations: 4

plot_model(cost75.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

cost100.B.drift.outcomes <- rbind(drift5000[[2]], s1.new.B0.cost100[[2]], s1.new.B.1.cost100[[2]], s1.new.B.2.cost100[[2]], s1.new.B.3.cost100[[2]], s1.new.B.4.cost100[[2]], s1.new.B.5.cost100[[2]], s1.new.B.6.cost100[[2]], s1.new.B.7.cost100[[2]], s1.new.B.8.cost100[[2]], s1.new.B.9.cost100[[2]], s1.new.B1.cost100[[2]])

cost100.B.drift.outcomes$B <- NA
cost100.B.drift.outcomes$B[1:5000] <- "drift"
cost100.B.drift.outcomes$B[5001:10500] <- rep(seq(0,1,0.1), each=500)
cost100.B.drift.outcomes$B <- factor(cost100.B.drift.outcomes$B, levels=c("drift","0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1"))

library(colorRamps)
cost100.B.drift.outcomes$green.outcome <- NA
cost100.B.drift.outcomes$green.outcome[cost100.B.drift.outcomes$e.var %in% blue2green(10)[1:5]] <- FALSE
cost100.B.drift.outcomes$green.outcome[cost100.B.drift.outcomes$e.var %in% blue2green(10)[6:10]] <- TRUE

library(sjPlot)
cost100.model <- glm(green.outcome~B, data=cost100.B.drift.outcomes, family="binomial")
summary(cost100.model)

Call:
glm(formula = green.outcome ~ B, family = "binomial", data = cost100.B.drift.outcomes)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5877  -1.1603   0.8166   1.1946   1.4336  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04040    0.02886  -1.400  0.16165    
B0           0.94810    0.11326   8.371  < 2e-16 ***
B0.1         0.96734    0.11531   8.389  < 2e-16 ***
B0.2         0.76707    0.10987   6.981 2.92e-12 ***
B0.3         0.57884    0.10673   5.423 5.84e-08 ***
B0.4         0.45007    0.10660   4.222 2.42e-05 ***
B0.5         0.37937    0.10415   3.642  0.00027 ***
B0.6         0.11600    0.10451   1.110  0.26704    
B0.7        -0.14668    0.10624  -1.381  0.16740    
B0.8        -0.33584    0.10603  -3.167  0.00154 ** 
B0.9        -0.27913    0.10483  -2.663  0.00775 ** 
B1          -0.54423    0.10899  -4.994 5.93e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12743  on 9202  degrees of freedom
Residual deviance: 12437  on 9191  degrees of freedom
  (1297 observations deleted due to missingness)
AIC: 12461

Number of Fisher Scoring iterations: 4

plot_model(cost100.model, show.values = TRUE, value.offset=.3, title="Odds ratios for green outcome")

#model 1 cost model outputs

B.model1.cost0.outputs <- data.frame(model = rep(NA, 11), cost = rep(0,11), B = seq(0,1,0.1), OR = exp(cost0.model$coefficients)[-1], low = exp(confint(cost0.model))[-1,1], high = exp(confint(cost0.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.cost10.outputs <- data.frame(model = rep(NA, 11), cost = rep(10,11), B = seq(0,1,0.1), OR = exp(cost10.model$coefficients)[-1], low = exp(confint(cost10.model))[-1,1], high = exp(confint(cost10.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.cost25.outputs <- data.frame(model = rep(NA, 11), cost = rep(25,11), B = seq(0,1,0.1), OR = exp(cost25.model$coefficients)[-1], low = exp(confint(cost25.model))[-1,1], high = exp(confint(cost25.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.cost50.outputs <- data.frame(model = rep(NA, 11), cost = rep(50,11), B = seq(0,1,0.1), OR = exp(cost50.model$coefficients)[-1], low = exp(confint(cost50.model))[-1,1], high = exp(confint(cost50.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.cost75.outputs <- data.frame(model = rep(NA, 11), cost = rep(75,11), B = seq(0,1,0.1), OR = exp(cost75.model$coefficients)[-1], low = exp(confint(cost75.model))[-1,1], high = exp(confint(cost75.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

B.model1.cost100.outputs <- data.frame(model = rep(NA, 11), cost = rep(100,11), B = seq(0,1,0.1), OR = exp(cost100.model$coefficients)[-1], low = exp(confint(cost100.model))[-1,1], high = exp(confint(cost100.model))[-1,2])

Waiting for profiling to be done...
Waiting for profiling to be done...

cost.model1.outputs <- rbind(B.model1.cost0.outputs, B.model1.cost10.outputs, B.model1.cost25.outputs, B.model1.cost50.outputs, B.model1.cost75.outputs, B.model1.cost100.outputs)

cost.model1.outputs$model[1:11] <- "cost = 0"
cost.model1.outputs$model[12:22] <- "cost = 10"
cost.model1.outputs$model[23:33] <- "cost = 25"
cost.model1.outputs$model[34:44] <- "cost = 50"
cost.model1.outputs$model[45:55] <- "cost = 75"
cost.model1.outputs$model[56:66] <- "cost = 100"


cost.model1.outputs$model <- factor(cost.model1.outputs$model)

ggplot(data = cost.model1.outputs, aes(x = as.numeric(as.character(B)), y = OR, color = model))+
geom_line()+
geom_point(size = 1)+
geom_ribbon(aes(y = OR, ymin = low, ymax = high, fill=model), alpha = 0.2)+
scale_color_manual(values=c(blue2green(6)))+
geom_hline(yintercept=1, linetype="dashed", color = "red")+
ylim(0, 4)+
ylab("Odds ratios for a green outcome")+
xlab("Value of B")+
scale_x_continuous(breaks = seq(0, 1, .1))+
theme(legend.title=element_blank())